#include "AnaDef.h"

const Double_t EPSILON = 1e-10;

const Int_t axJpt = 1;
const Int_t axPDG = 7;
Int_t axVar = -999;
Int_t fNrbX = -999;

void GetHist(const TList *ll, THnSparseD * &hh, const TString name)
{
  hh=(THnSparseD*)ll->FindObject(name);
  if(!hh){
    printf("no hh %s!!\n",name.Data()); exit(1);
  }
  printf("\nHh %s content:\n", name.Data());
  hh->Print();
  const Int_t ndim = hh->GetNdimensions();
  for(Int_t ii=0; ii<ndim; ii++){
    const TAxis *ax = hh->GetAxis(ii);
    const Int_t nbins = ax->GetNbins();
    printf("%s axis %d: %s, nbin= %d\n", name.Data(), ii, ax->GetTitle(), nbins);
  }
  printf("\n");
}

void OutputHist(THnSparseD * hh, const TString tag)
{
  const TString dirname(Form("histRec%s_LHC%s", hh->GetName(), tag.Data()));
  gSystem->Exec(Form("mkdir %s", dirname.Data()));
  TFile * fout = new TFile((dirname+"/")+(dirname+".root"),"recreate");

  TH2D * hpdg = hh->Projection(axPDG, axVar);
  hpdg->SetName("mcpdgVsPt");
  printf("test hpdg nx %d ny %d entries %e integral %e\n", hpdg->GetNbinsX(), hpdg->GetNbinsY(), hpdg->GetEntries(), hpdg->Integral(1, hpdg->GetNbinsX(), 1, hpdg->GetNbinsY()));
  if(fNrbX>1){
    printf("rebinning...\n");
    hpdg->RebinX(fNrbX);
  }

  TH1D *htype[]={0x0,0x0,0x0,0x0};
  for(Int_t itype=0; itype<4; itype++){
    htype[itype] = hpdg->ProjectionX(Form("mcpar%dPt",itype), itype+1, itype+1);
    htype[itype]->Write();
  }

  hpdg->Write();
  fout->Save();
  fout->Close();
  delete fout;
  
  //protected by TH1::AddDirectory(kFALSE);
  delete hpdg;
  for(Int_t ii=0; ii<4; ii++){
    delete htype[ii];
  }
}

void getOutputHist(THnSparseD * hh, const TString tag, const Bool_t kloopJpt)
{
  const Int_t nJpt = kloopJpt?  hh->GetAxis(axJpt)->GetNbins() : 1; 

  //remove over-/underflow of other bins
  SetAxisBit(hh, kTRUE);  
 
  TH1D * hall = hh->Projection(axJpt);
  for(Int_t ijpt=0; ijpt<=nJpt+1; ijpt++){
    const Double_t icont = hall->GetBinContent(ijpt);
    printf("%s Entries jpt %d: %e\n", hh->GetName(), ijpt, icont);
  }
  TFile * fall = new TFile(Form("%s_LHC%s_allEta_allJpt.root", hh->GetName(), tag.Data()),"recreate");
  hall->Write();
  fall->Save();
  fall->Close();
  delete fall;
  delete hall;

  //important, otherwise the 2D projection in OutputHist will have projects nothing.
  SetAxisBit(hh, kFALSE);

  if(kloopJpt){
    hh->GetAxis(axJpt)->SetRange(1, nJpt);
    //go after SetRange to deselect over-/under-flow
    //http://root.cern.ch/root/html528/THnSparse.html#THnSparse:IsInRange
    hh->GetAxis(axJpt)->SetBit(TAxis::kAxisRange);
  }

  OutputHist(hh, Form("%s_ETAall_JPTall",tag.Data()));

  if(!kloopJpt)
    return;
  
  for(Int_t ijpt = 1; ijpt<=nJpt; ijpt++){
    printf("Projecting jpt %d ...\n", ijpt);
    hh->GetAxis(axJpt)->SetRange(ijpt, ijpt);
    
    OutputHist(hh, Form("%s_ETAall_JPT%d", tag.Data(), ijpt-1));
  }

  printf("\n");

}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=5){
    printf("bad argc! %d\n", argc); 
    return 1;
  }
  const TString tag(argv[2]);
  TString var(argv[3]);
  var.ToUpper();
  const Int_t tpccutmode=atoi(argv[4]);

  const Int_t axPt  = 4;
  const Int_t axZ = 5;
  const Int_t axXi = 6;
  
  if(var=="PT"){
    axVar = axPt;
    fNrbX = 1;
  }
  else if(var=="Z"){
    axVar = axZ;
    fNrbX = 1;
  }
  else if(var=="XI"){
    axVar = axXi;
    fNrbX = 1;
  }
  else{
    printf("bad var! %s\n", var.Data()); return 1;
  }

  /* http://root.cern.ch/root/html/TH1.html
     When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
     h->SetDirectory(0);          for the current histogram h   
     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference  
     When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are auto \
     matically deleted.    */
  //very important
  TH1::AddDirectory(kFALSE);

  TFile *_file0 = new TFile(argv[1]);
  printf("\nFile content:\n");
  _file0->ls();

  TString dns[2];
  TString lns[2];

  dns[0]="PWGJE_IDFFTCF_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc0_cl0";
  dns[1]="PWGJE_IDFFTCF_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc1_cl0";
  
  lns[0]="idfftcf_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc0_cl0";
  lns[1]="idfftcf_clustersAOD_ANTIKT04_B0_Filter00272_Cut00150_Skip00_clustersAODMC2_ANTIKT04_B0_Filter00272_Cut00150_Skip00_AODMC2b_AODMCb_tpc1_cl0";

  TDirectory *fdir=_file0->GetDirectory(dns[tpccutmode]);
  if(!fdir){
    printf("no fdir!!\n");exit(1);
  }
  printf("\nDirectory content:\n");
  fdir->ls();

  TList * ll=(TList*)fdir->Get(lns[tpccutmode]);
  if(!ll){
    printf("no list!!\n");exit(1);
  }
  printf("\nList content:\n");
  ll->ls();

  THnSparseD * hhITPC=0x0;
  THnSparseD * hhIdff=0x0;

  printf("\nGetHist ...\n");
  GetHist(ll, hhIdff, "IDFFHist");
  GetHist(ll, hhITPC, "InclusiveTPC");

  printf("\ngetOutputHist ...\n");
  getOutputHist(hhIdff, tag, 1);
  getOutputHist(hhITPC, tag, 0);

  return 0;
}
